\documentclass{article}

\usepackage[top=25mm,bottom=25mm,left=25mm,right=25mm]{geometry}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{subfig}
\usepackage{epstopdf}
\usepackage[framed,autolinebreaks,useliterate]{mcode}
\setlength{\parindent}{0pt}

\begin{document}

\title{ESP - Digital lab 4}
\author{Floris van Nee \& Simon Dirlik}

\maketitle
\section{} %4.1
The desired values are calculated with following matlab code:
\begin{lstlisting}
	F =   [0.1 0.18 0.62 0.70];
	A =   [0 1 0];
	DEV = [0.001 1 0.001]; % 10^(0.2/20)/10^(60/20) = 0.001
	[N,Wn,BTA,FILTYPE] = kaiserord(F,A,DEV);
\end{lstlisting}
This gives $N=91$, however 91 was not the minimum order to get the desired transfer function. After trial and error we found that $N=85$ was sufficient. These values can be used in the function \verb|fir1|, which gives the transfer function (all b-coefficients in case of a FIR filter). The impulse response and the frequency response can be calculated and plotted with the function \verb|freqz| and \verb|impz|. The results of the following matlab code are shown in figures \ref{fig:impz41} and \ref{fig:freqz41}.
\begin{lstlisting}
	B = fir1(N-6, Wn, FILTYPE, kaiser( N+1-6,BTA ),'noscale'); % N-6=85;
	figure;	impz(B,1); % impuls response
	figure;	freqz(B); % frequency response
\end{lstlisting}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{impz41.eps}
	\end{center}
	\caption{The impulse response computed using fir1 and impz.}
	\label{fig:impz41}
\end{figure}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{freqz41.eps}
	\end{center}
	\caption{The frequency response computed using fir1 and freqz.}
	\label{fig:freqz41}
\end{figure}

\section{}%4.2
Following the same methods as in 4.1, but using \verb|firpm| and \verb|firpmord| in stead of \verb|fir1| and \verb|kaiserord|, we get the following matlab code, the results of this code are shown in figure \ref{fig:impz42} and \ref{fig:freqz42}.\\Note: In 4.1 the estimator function gave a value for the order which was higher than the minimal correct value, in this exercise the value for \verb|n| was to low.
\begin{lstlisting}
	[n,fo,mo,w] = firpmord(F,A,DEV);
	b = firpm(n+7,fo,mo,w); % n=25; not enough
	figure;	impz(b,1); % impuls response
	figure;	freqz(b); % frequency response
\end{lstlisting}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{impz42.eps}
	\end{center}
	\caption{The impulse response computed using firpm and impz.}
	\label{fig:impz42}
\end{figure}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{freqz42.eps}
	\end{center}
	\caption{The frequency response computed using firpm and freqz.}
	\label{fig:freqz42}
\end{figure}

\section{}%4.3
The FIR filter using the Kaiser window gives a function of higher order than the optimal equiripple function. This results in a longer impulse response and smoother lines in the magnitude plot, which is expected when a higher order is used. The advantage of the equiripple is that it requires less calculations while still staying within the boundaries that were specified for the ripple.
\clearpage
\section{}%4.4
	The function \verb|bitround(b,bit)| rounds the elements of vector \verb|b| to $2^{bit}$ levels, so it can be used to quantize a signal. The results of the following matlab code are shown in figures \ref{fig:freqz44b16}, \ref{fig:freqz44b12} and \ref{fig:freqz44b10}.
\begin{lstlisting}
	bits=10; % we did this for bits = 10, 12 and 16
	figure;	freqz(bitround(b,bits)); % frequency response
\end{lstlisting}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{freqz44b16.eps}
	\end{center}
	\caption{The frequency response of the equiripple filter rounded to 16 bits.}
	\label{fig:freqz44b16}
\end{figure}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{freqz44b12.eps}
	\end{center}
	\caption{The frequency response of the equiripple filter rounded to 12 bits.}
	\label{fig:freqz44b12}
\end{figure}
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{freqz44b10.eps}
	\end{center}
	\caption{The frequency response of the equiripple filter rounded to 10 bits.}
	\label{fig:freqz44b10}
\end{figure}
Rounding of causes error in the frequency response. The figures show that 16 bits accuracy is enough create a filter which meets the specification, with 10 and 12 bit accuracy the specification is no longer met. The stopband parts are now less than 60dB from the passband part. We can conclude from this that the digital filter needs a high accuracy for the calculation to be valid.
\clearpage
\section{}%4.5
To cancel out H, we want a frequency response of H inverse:
\[\begin{split}
	H^{-1}	& = \dfrac {1} {\dfrac {sin(\omega\dfrac{T_s}{2})} {\omega\dfrac{T_s}{2}}} \\
			& = \dfrac {\omega\dfrac{T_s}{2}} {sin(\omega\dfrac{T_s}{2})}\\
\end{split}\]
The IDFT can be calculated using the matlab function \verb|ifft|. We can compare the frequency response of the result of that function with the original desired frequency response. This is done with the following matlab code.
\begin{lstlisting}
	order=20; n=(2*order)+1; step=2*pi/(n-1);
	omts = -pi:step:pi;
	Hi = (omts./2)./(sin(omts./2)); % the desired frequency response H inverse
	Hi(order+1)=1; % to get rid of the NaN at the zero-division.
	idft = ifft(Hi);
	computedfr=freqz(idft,1,n);
	difference = abs(computedfr')-abs(Hi);
	figure; plot(abs(computedfr));
	figure; plot(abs(Hi));
	figure; plot(difference);
\end{lstlisting}
The results are shown in figures \ref{fig:denc45} and \ref{fig:diff45}. The figures show that the computed frequency is not completely accurate.
% wat kan ik hier nou over zeggen.
\begin{figure}[h]
	\begin{center}
		\subfloat[The desired frequency response.]{\label{fig:desired45}\includegraphics[width=0.5\textwidth]{desired45.eps}}
		\subfloat[The computed frequency response.]{\label{fig:computed45}\includegraphics[width=0.5\textwidth]{computed45.eps}}
	\end{center}
	\caption{}
	\label{fig:denc45}
\end{figure}
\clearpage
\begin{figure}[h]
	\begin{center}
		\includegraphics[width=\textwidth]{diff45.eps}
	\end{center}
	\caption{The difference between the frequency responses.}
	\label{fig:diff45}
\end{figure}
\section{}%4.6
Matlab has functions to estimate the order for all of these filters, the following code shows how to use them:
\begin{lstlisting}
	% specifications
	Wp = [0.18 0.62]; Ws = [0.1 0.7];
	Rp = 0.2; Rs = 60;

	% Butterworth filter
	[n,Wn] = buttord(Wp,Ws,Rp,Rs); % n=22
	[b,a] = butter(n,Wn);

	% Chebyshev (type 1) filter
	[n,Wn] = cheb1ord(Wp,Ws,Rp,Rs); % n=10
	[b,a] = cheby1(n,Rp,Wn);

	% Cauer (elliptic) filter
	[n,Wn] = ellipord(Wp,Ws,Rp,Rs); % n=7
	[b,a] = ellip(n,Rp,Rs,Wn); % transfer function sum(b(i)/1+a(i))
\end{lstlisting}
In a previous exercise the order gotten from the estimation functions were not the correct minimal orders, in this case we found that the functions do give the correct orders.
\section{}%4.7
There are matlab functions to calculate and plot the impulse response, step response and the groupdelay. The poles and zeros can be calculated in 2 ways, either way the result is a vector of complex numbers. To plot the poles and zeros in the z-plane we need to transform to polar form. This is shown in following the code:
\begin{lstlisting}
	% specifications
	Wp = [0.18 0.62]; Ws = [0.1 0.7];
	Rp = 0.2; Rs = 60;

	% Cauer (elliptic) filter
	[n,Wn] = ellipord(Wp,Ws,Rp,Rs); % n=7
	[b,a] = ellip(n,Rp,Rs,Wn); % transfer function sum(b(i)/1+a(i))
	figure;	impz(b,a); % impulse response
	figure;	stepz(b,a); % step response
	figure;	grpdelay(b,a); % matlab group delay function
	[z,p,k] = ellip(n,Rp,Rs,Wn); % poles and zeros, k is gain
	% or: z = roots(b); p = roots(a);
	figure; grid ON;
	polar(atan(imag(z)./real(z)) , abs(z), 'o'); % tan(theta) = opposite/adjacent
	hold;
	polar(atan(imag(p)./real(p)) , abs(p), 'x'); % theta = atan(imaginary/real)
	hold;
\end{lstlisting}
The results are shown in figure \ref{fig:all47}. The pole-zero diagram show that all the zeros are on the unit circle, which is expected because the unit circle is the imaginary axis in the s-plane. It also show that all the poles are within the unit circle which means the system is stable, because that part is the left hand side of the s-plane. The minimum orders are: Butterworth:22, Chebyshev:10 and Cauer:7.
\begin{figure}[h]
	\begin{center}
		\subfloat[The frequency response.]{\label{fig:freqz47}\includegraphics[width=0.5\textwidth]{freqz47.eps}}
		\subfloat[The group delay time.]{\label{fig:grpdelay47}\includegraphics[width=0.5\textwidth]{grpdelay47.eps}}\\
		\subfloat[The pole-zero diagram. (0=zero, x=pole)]{\label{fig:penz47}\includegraphics[width=0.5\textwidth]{penz47.eps}}
		\subfloat[The impulse response.]{\label{fig:impz47}\includegraphics[width=0.5\textwidth]{impz47.eps}}\\
		\subfloat[The step response.]{\label{fig:step47}\includegraphics[width=0.5\textwidth]{step47.eps}}
	\end{center}
	\caption{Cauer filter.}
	\label{fig:all47}
\end{figure}
\section{}%4.8
Rounding can be done with \verb|bitround|, the rest is similar to exercise 4.7
\begin{lstlisting}
	br=bitround(b,12);
	ar=bitround(a,12);
	figure;	freqz(br,ar);
	zr = roots(b);
	pr = roots(a);
	figure;	grid ON;
	polar(atan(imag(zr)./real(zr)) , abs(zr), 'o'); % tan(theta) = overstaand/aanliggend
	hold;
	polar(atan(imag(pr)./real(pr)) , abs(pr), 'x'); % theta = atan(imaginair/reeel)
	hold;
\end{lstlisting}
The results are shown in figure \ref{fig:all48}. It shows that the frequency response is affected like in exercise 4.4, but the poles and zeros are exactly the same.
\begin{figure}[h]
	\begin{center}
		\subfloat[The frequency response.]{\label{fig:freqz48}\includegraphics[width=0.5\textwidth]{freqz48.eps}}
		\subfloat[The pole-zero diagram. (0=zero, x=pole)]{\label{fig:penz48}\includegraphics[width=0.5\textwidth]{penz48.eps}}
	\end{center}
	\caption{Cauer filter rounded off.}
	\label{fig:all48}
\end{figure}
\section{}%4.9
The orders of the IIR filters are generally lower than the orders of the FIR filters, This does give the FIR filters smoother lines. Rounding is hard for both types of filters.
% wat nog ...
\end{document}